function [ Ay ] = getAy( i,j,U,V,NJve,NIve ,NIvt,h)

V_cen  = V ( Vmap (i  ,j  ) );  % center of finite volume
if (i==NIvt)   V_right = -V_cen; else V_right = V ( Vmap (i+1,j ) ); end

if (j<=NJve)
    % Right and left boundaries, considering ghost nodes
    if (i==1) V_left = -V_cen; else V_left = V ( Vmap (i-1,j ) ); end
else
    % Right and left boundaries, considering ghost nodes
    if (i==NIve+1) V_left  = -V_cen; else V_left  = V ( Vmap (i-1,j ) ); end
end
V_bot = V ( Vmap (i,j-1));
V_top = V ( Vmap (i,j+1));

U_tr    = U ( Umap ( i+1,j  ) );
U_tl    = U ( Umap ( i  ,j  ) );
U_br    = U ( Umap ( i+1,j-1) );
U_bl    = U ( Umap ( i  ,j-1) );

Ay  = + ( U_tr + U_br ) * ( V_cen + V_right ) ...
    + ( V_cen + V_top )^2 ...
    - ( U_tl + U_bl ) * ( V_cen + V_left ) ...
    - ( V_cen + V_bot )^2;

Ay  = Ay/(4*h);

end

